import numpy as np
import matplotlib.pyplot as plt
from pylab import imshow,gray,show
from mpl_toolkits import mplot3d
import re
import matplotlib.cm as cm
from matplotlib.colors import Normalize
import matplotlib.pyplot as plotLine
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import LinearRegression
import statsmodels.api as sm
import statistics as stat
from sklearn.metrics import mean_squared_error
import itertools
import math
import matplotlib.transforms as mtransforms
import matplotlib.pyplot as plt
from scipy.signal import savgol_filter
import os
from sklearn import preprocessing
def get_time_dipole(name, dcdfreq, flaginittime):
data=np.loadtxt(f'{name}',dtype=float)
#print(data)
time=[]
dipolez=[]
anglex=[]
angley=[]
anglez=[]
for row in data:
if flaginittime==1:
time.append((row[0]*dcdfreq*2e-6)+0.005)
else:
time.append(row[0]*dcdfreq*2e-6)
#time.append(row[0])
#dipole.append(row[4])
dipolez.append(row[3])
anglex.append(math.degrees(math.acos(row[1]/row[4])))
angley.append(math.degrees(math.acos(row[2]/row[4])))
anglez.append(math.degrees(math.acos(row[3]/row[4])))
mean_anglez=stat.mean(anglez)
st_dev=stat.stdev(anglez)
#print(name,"mean:", mean_anglez,"st_dev:",st_dev)
return time, dipolez, anglex, angley, anglez, mean_anglez , st_dev
def get_average_anglez(data,init,finish,name):
data_sel=[]
#selet data
for i in range(init,finish+1):
data_sel.append(data[i])
mean_anglez=stat.mean(data_sel)
st_dev=stat.stdev(data_sel)/(len(mean_anglez))**0.5
print(name,"mean:", mean_anglez,"st_dev:",st_dev)
return mean_anglez , st_dev
def get_average_data(Pz,theta,init,finish,name):
Pz_sel=[]
theta_sel=[]
#selet data
for i in range(init,finish):
Pz_sel.append(Pz[i])
theta_sel.append(theta[i])
mean_Pz=stat.mean(Pz_sel)
st_dev_Pz=stat.stdev(Pz_sel)
mean_theta=stat.mean(theta_sel)
st_dev_theta=stat.stdev(theta_sel)
print(name)
print("mean_Pz:",mean_Pz,"st_dev_Pz:",st_dev_Pz)
print("mean_theta:",mean_theta,"st_dev_theta:",st_dev_theta)
num_bins=20
# plt.hist(theta, num_bins)
# plt.show()
# plt.hist(Pz, num_bins)
# plt.show()
return mean_Pz, st_dev_Pz,mean_theta,st_dev_theta
def write_point(x, y, ax,hpos, vpos, xpos, ypos):
x=round(x, 3)
y=round(y, 3)
trans_offset = mtransforms.offset_copy(ax.transData, fig=ax, x=xpos, y=ypos, units='dots')
ax.text(x, y, f' {x}ns', verticalalignment=f'{vpos}', horizontalalignment=f'{hpos}',fontsize=18,zorder=2,transform=trans_offset) #put fontsize 11 por DW-0.2V/nm and another case 18
ax.scatter(x, y, color='black', marker='o',zorder=4)
def get_time_energy(name, dcdfreq):
data=np.loadtxt(f'{name}',dtype=float)
#print(data)
time=[]
energy=[]
for row in data:
#time.append(row[0]*dcdfreq*2e-6)
time.append((row[0]*dcdfreq*2e-6)+0.005)
#time.append(row[0])
energy.append(row[1])
return time, energy
def compute_max_normalzD8C(field,name_sim):
maxsPw=[]
maxsPm=[]
sim=[f'{location}{name_sim}_ZDMW.dat',f'{location}{name_sim}.dat']
for i in field:
sim.append (f'{location}{name_sim}0{i}E_ZDMW.dat')
sim.append(f'{location}{name_sim}0{i}E.dat')
for enum,n in enumerate(sim):
if enum%2==0:
maxsPw.append(max(get_time_dipole(n,2500,0)[1]))
else:
maxsPm.append(max(get_time_dipole(n,2500,0)[1]))
Pmax=[max(maxsPw),max(maxsPm)]
print("Pz max:", Pmax)
return(Pmax)
Problem:
What parameter allows monitoring the membrane permeabilization state and how is it used?
Answer:
The key parameter for monitoring the membrane permeabilization state is
the dipole moment. By observing the behavior of the z-axis component of
the membrane dipole moment, , and the angle between the membrane dipole vector and the z-axis, , we observe a drastic change as membrane permeabilization occurs.
Additionally, we identify the time at which the first molecules of water
enter the interior of the membrane, marking the beginning of
permeabilization (). These values were plotted with respect to the magnitude of the electric field.
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
Fntsize=18
#plot settings
fig = plt.figure()
fig.set_figheight(14)
fig.set_figwidth(20)
ax = plt.subplot2grid(shape=(2, 2), loc=(0, 0), colspan=2)
ax1 = plt.subplot2grid(shape=(2, 2), loc=(1, 0), colspan=1)
ax2 = plt.subplot2grid(shape=(2, 2), loc=(1, 1), colspan=1)
colors = cm.rainbow(np.linspace(0, 1, 7))
colors = np.arange(0, 1, 5)
colors=cm.viridis(np.linspace(0, 1, 7))
#files Z-DMw
DMPC=f'{location}DMPC.dat' #0V/nm
DMPCE2=f'{location}2DMPCE.dat' #0.1V/nm
DMPCE0=f'{location}0DMPCE.dat' #0.2V/nm
#DMPCE0c=f'{location}0DMPCE.dat' #0.2V/nm
DMPCE4=f'{location}4DMPCE.dat' #0.3V/nm
DMPCE5=f'{location}5DMPCE.dat' #0.4V/nm
DMPCE1=f'{location}1DMPCE.dat' #0.5V/nm
#Z-DMw
ax.plot(get_time_dipole(DMPC,50000,1)[0],get_time_dipole(DMPC,50000,1)[1],label='DW-0V/nm',color=colors[0])
ax.plot(get_time_dipole(DMPCE2,2500,1)[0],get_time_dipole(DMPCE2,2500,1)[1],label='DW-0.1V/nm',color=colors[1])
#ax.plot(get_time_dipole(DMPCE3,2500,1)[0],get_time_dipole(DMPCE3,2500,1)[1],label='DW-0.16V/nm',color=colors[2])
ax.plot(get_time_dipole(DMPCE0,2500,1)[0],get_time_dipole(DMPCE0,2500,1)[1],label='DW-0.2V/nm',color=colors[2])
#ax.plot(get_time_dipole(DMPCE0c,2500,1)[0],get_time_dipole(DMPCE0c,2500,1)[1],label='DWc-0.2V/nm',color='black')
ax.plot(get_time_dipole(DMPCE4,2500,1)[0],get_time_dipole(DMPCE4,2500,1)[1],label='DW-0.3V/nm',color=colors[3])
ax.plot(get_time_dipole(DMPCE5,2500,1)[0],get_time_dipole(DMPCE5,2500,1)[1],label='DW-0.4V/nm',color=colors[4])
ax.plot(get_time_dipole(DMPCE1,2500,1)[0],get_time_dipole(DMPCE1,2500,1)[1],label='DW-0.5V/nm',color=colors[5])
#angle of the dipole moment with respect to the z plane.
ax1.plot(get_time_dipole(DMPC,50000,0)[0],get_time_dipole(DMPC,50000,1)[4],label='DW-0V/nm',color=colors[0])
ax1.plot(get_time_dipole(DMPCE2,2500,1)[0],get_time_dipole(DMPCE2,2500,1)[4],label='DW-0.1V/nm',color=colors[1])
#ax1.plot(get_time_dipole(DMPCE3,2500,1)[0],get_time_dipole(DMPCE3,2500,1)[4],label='DW-0.16V/nm',color=colors[2])
ax1.plot(get_time_dipole(DMPCE0,2500,1)[0],get_time_dipole(DMPCE0,2500,1)[4],label='DW-0.2V/nm',color=colors[2])
#x1.plot(get_time_dipole(DMPCE0c,2500,1)[0],get_time_dipole(DMPCE0c,2500,1)[4],label='DWc-0.2V/nm',color='black')
ax1.plot(get_time_dipole(DMPCE4,2500,1)[0],get_time_dipole(DMPCE4,2500,1)[4],label='DW-0.3V/nm',color=colors[3])
ax1.plot(get_time_dipole(DMPCE5,2500,1)[0],get_time_dipole(DMPCE5,2500,1)[4],label='DW-0.4V/nm',color=colors[4])
ax1.plot(get_time_dipole(DMPCE1,2500,1)[0],get_time_dipole(DMPCE1,2500,1)[4],label='DW-0.5V/nm',color=colors[5])
#plot names and legend
ax.set_xlabel('Time(ns)',fontsize=Fntsize)
ax.set_ylabel('$P_{zDMPC}$ (D)',fontsize=Fntsize)
ax1.set_ylabel('$\Theta_{zDMPC}$ (Degree)',fontsize=Fntsize)
ax1.set_xlabel('Time(ns)',fontsize=Fntsize)
ax.set_title('A)',fontsize=Fntsize)
ax1.set_title('B)',fontsize=Fntsize)
ax1.legend(loc="upper right",fontsize=Fntsize,ncol=2)
#ax.legend(loc="center right")
plt.tight_layout(pad=2)
ax.set_xlim(0,4)
ax1.set_xlim(0,4)
tp=[0,191,36,25,5]
#Put points
#0.2V/nm
#write_point((get_time_dipole(DMPCE0,2500,1)[0])[145],(get_time_dipole(DMPCE0,2500,1)[1])[125],ax,'right')
write_point((get_time_dipole(DMPCE0,2500,1)[0])[191],(get_time_dipole(DMPCE0,2500,1)[1])[191],ax,'left','center',-15,13)
#0.3V/nm
write_point((get_time_dipole(DMPCE4,2500,1)[0])[36],(get_time_dipole(DMPCE4,2500,1)[1])[36],ax,'left','center',5,0)
#0.4V/nm
write_point((get_time_dipole(DMPCE5,2500,1)[0])[25 ],(get_time_dipole(DMPCE5,2500,1)[1])[25],ax,'left','center',-4,-10)
#0.5 v/nm
write_point((get_time_dipole(DMPCE1,2500,1)[0])[5],(get_time_dipole(DMPCE1,2500,1)[1])[5],ax,'left', 'center',-8, 24)
tp=[(get_time_dipole(DMPCE0,2500,1)[0])[191],(get_time_dipole(DMPCE4,2500,1)[0])[36],(get_time_dipole(DMPCE5,2500,1)[0])[25],(get_time_dipole(DMPCE1,2500,1)[0])[5]]
E=[0.2,0.3,0.4,0.5]
ax2.scatter(E,tp)
ax2.plot(E,tp, color='black')
ax2.set_title('C)',fontsize=Fntsize)
ax2.set_ylabel('$t_{p} (ns)$',fontsize=Fntsize)
ax2.set_xlabel('$E$ $(V/nm)$',fontsize=Fntsize)
plt.tight_layout()
plt.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Z-DMW_DW.png')
print((tp[0]-tp[1]))
print((tp[1]-tp[2]))
0.7749999999999999 0.05499999999999999
Def variables avg dipole
AvgDM_DMPC_0={}
AvgDM_Water_0={}
AvgDM_DMPC_01={}
AvgDM_Water_01={}
# fig, axs = plt.subplots(1, 1,figsize=(16, 8))
DMPCE0=f'{location}0DMPCE_ZDMW.dat' #0.2V/nm
DMPCE0o=f'{location}0DMPCE_ZDMW_origin.dat' #0.2V/nm
DMPCE0e=f'{location}0DMPCE_ZDMW_e.dat' #0.2V/nm
# ax=axs
fig, axes = plt.subplots(3, 4, figsize=(10, 10))
fig.set_figheight(14)
fig.set_figwidth(15)
axes = axes.flatten()
axes[0] = plt.subplot2grid(shape=(3, 4), loc=(0, 0), colspan=3)
#Compute derivate of the dipole moment
xvalues=get_time_dipole(DMPCE0,2500,1)[0]
yvalues=get_time_dipole(DMPCE0,2500,1)[1]
#xvalues = savgol_filter(xvalues, 9, 2)
yvalues = savgol_filter(yvalues, 41, 4)
dZ_DMW=np.gradient(yvalues,xvalues)
axes[0].plot(xvalues,dZ_DMW,label=r'$\frac{dP_{zw}}{dt}$',color='gray')
axes[0].plot(get_time_dipole(DMPCE0,2500,1)[0],get_time_dipole(DMPCE0,2500,1)[1],label=r'$P_{zw}$',color=colors[2])
axes[0].set_xlabel('Time(ns)',fontsize=Fntsize)
axes[0].set_ylabel(r'$P_{zw} (D) \ , \ \frac{dP_{zw}}{dt} (D/ns)$',fontsize=Fntsize)
ax=axes[0]
ax.set_xlim(0,4)
#ax.set_ylim(-50000,50000)
ax.legend(loc="upper right",fontsize=Fntsize,ncol=2)
times=[0,122,191,241,268,304,332,441,799]
a=3
for t in times:
img=plt.imread(f'/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMPCE/0DMPCE/0DMPCE_{t}.png')
axes[a].imshow(img)
#axes[a].axis('off') # Turn off axis labels
tns=round((t*2500*2e-6)+0.005, 3)
axes[a].set_title(f'{tns}ns',fontsize=Fntsize)
axes[a].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[a].xaxis.set_visible(False)
#axes[4].yaxis.set_visible(False)
axes[3].set_ylabel('E',fontsize=Fntsize)
axes[a].tick_params(left=False, labelleft=False)
plt.setp(axes[a].spines.values(), visible=False)
a+=1
start_point = (0, 556)
end_point = (0, 100)
axes[3].arrow(start_point[0], start_point[1],
end_point[0] - start_point[0], end_point[1] - start_point[1],
color='black', width=3, head_width=10, head_length=10)
axes[3].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[3].xaxis.set_visible(False)
#axes[4].yaxis.set_visible(False)
axes[3].set_ylabel('$E_{z}=0.2V/nm$',fontsize=20)
axes[3].tick_params(left=False, labelleft=False)
plt.setp(axes[3].spines.values(), visible=False)
# ax.plot(get_time_dipole(DMPCE0o,2500,1)[0],get_time_dipole(DMPCE0o,2500,1)[1],label='DW-0.2V/nm_origin',color=colors[5])
# ax.plot(get_time_dipole(DMPCE0e,2500,1)[0],get_time_dipole(DMPCE0e,2500,1)[1],label='DW-0.2V/nm_e',color=colors[6])
#Put points
#0.2V/nm
write_point((get_time_dipole(DMPCE0,2500,1)[0])[0],(get_time_dipole(DMPCE0,2500,1)[1])[0],axes[0],'right','center',33,-7)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[122],(get_time_dipole(DMPCE0,2500,1)[1])[122],ax,'right','center',10,-13)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[191],(get_time_dipole(DMPCE0,2500,1)[1])[191],ax,'right','center',10,-13)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[241],(get_time_dipole(DMPCE0,2500,1)[1])[241],ax,'right','center',10,-10)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[268],(get_time_dipole(DMPCE0,2500,1)[1])[268],ax,'left','center',3,-4)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[304],(get_time_dipole(DMPCE0,2500,1)[1])[304],ax,'left','center',5,0)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[332],(get_time_dipole(DMPCE0,2500,1)[1])[332],ax,'left','center',5,0)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[441],(get_time_dipole(DMPCE0,2500,1)[1])[441],ax,'left','center',1,8)
write_point((get_time_dipole(DMPCE0,2500,1)[0])[799],(get_time_dipole(DMPCE0,2500,1)[1])[799],ax,'left','center',-30,-10)
#Adjust layout to prevent overlapping
plt.tight_layout()
plt.show()
show()
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Z-DMW_DW_0.2.png',dpi=300)
from PIL import Image
def conv_frame_ns(frames,dcdfreq):
timesns=[]
for f in frames:
timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
return(timesns)
loctp='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/tps/'
loctp20='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/tps20/'
loctpend='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/the end/'
loctpendtop='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/theEndtop/'
locations=[loctp,loctp20,loctpend,loctpendtop]
#psimage.save('myImage.png')
# fig, axs =pl
fig, axes = plt.subplots(4, 4, figsize=(10, 10))
fig.set_figheight(20)
fig.set_figwidth(20)
row=0
tpDMPCE0=[191,46,35,5]
simulations=[0,4,5,1]
#simulations=[0]
timetpns=conv_frame_ns(tpDMPCE0,2500)
timetp20ns=conv_frame_ns([fr+20 for fr in tpDMPCE0],2500)
timetpendns=conv_frame_ns([799,799,231,253],2500)
alltimens=[timetpns,timetp20ns,timetpendns,timetpendns]
fieldname=['0.2V/nm','0.3V/nm','0.4V/nm','0.5V/nm']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)
name='DMPCE'
for row,loc in enumerate(locations):
timesns=alltimens[row]
for col,idx in enumerate(simulations):
psimage=Image.open(f'{loc}{idx}{name}')
axes[row,col].imshow(psimage)
axes[row,col].set_title(f'{timesns[col]}ns',fontsize=Fntsize+5)
axes[row,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[row,col].xaxis.set_visible(False)
axes[row,col].tick_params(left=False, labelleft=False)
plt.setp(axes[row,col].spines.values(), visible=False)
#axes[row,col].yaxis.set_visible(False)
#plot field vector
if row==0:
axes[row,col].arrow(start_point[0], start_point[1],
end_point[0] - start_point[0], end_point[1] - start_point[1],
color='black', width=6, head_width=15, head_length=10)
axes[row,col].set_ylabel(f'$E_{z}$={fieldname[col]}',fontsize=20)
fig.suptitle('DW',fontsize=30)
# #plt.subplots_adjust(left=0.1,
# bottom=0.1,
# right=0.9,
# top=0.9,
# wspace=0.1,
# hspace=0)
plt.tight_layout(pad=0.5)
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/timetp_DW.png',dpi=300)
from PIL import Image
def conv_frame_ns(frames,dcdfreq):
timesns=[]
for f in frames:
timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
return(timesns)
loctin='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMin/'
loctp='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMtps/'
loctp20='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMtps20/'
loctpend='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMend/'
loctpendtop='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/theEndtop/'
locations=[loctin, loctp,loctp20,loctpend]
#psimage.save('myImage.png')
# fig, axs =pl
fig, axes = plt.subplots(2, 4, figsize=(10, 10))
fig.set_figheight(10)
fig.set_figwidth(15)
row=0
tp=[191,5]
simulations=[0,1]
#simulations=[0]
timetpns=conv_frame_ns(tp,2500)
timetp20ns=conv_frame_ns([fr+20 for fr in tp],2500)
timetpendns=conv_frame_ns([799,253],2500)
alltimens=[[0.005]*8,timetpns,timetp20ns,timetpendns,timetpendns]
fieldname=['0.2V/nm','0.5V/nm']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)
name='DMPCE'
for row,loc in enumerate(locations):
timesns=alltimens[row]
for col,idx in enumerate(simulations):
psimage=Image.open(f'{loc}{idx}{name}')
axes[col,row].imshow(psimage)
axes[col,row].set_title(f'{timesns[col]}ns',fontsize=Fntsize+5)
axes[col,row].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[col,row].xaxis.set_visible(False)
axes[col,row].tick_params(left=False, labelleft=False)
plt.setp(axes[col,row].spines.values(), visible=False)
if row==0:
axes[col,row].arrow(start_point[0], start_point[1],
end_point[0] - start_point[0], end_point[1] - start_point[1],
color='black', width=6, head_width=20, head_length=10)
axes[col,row].set_ylabel(f'$E_{z}$={fieldname[col]}',fontsize=20)
fig.suptitle('DW',fontsize=30)
plt.tight_layout(pad=0.5)
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Pdmpcarrowtimetp.png',dpi=300)
from PIL import Image
def conv_frame_ns(frames,dcdfreq):
timesns=[]
for f in frames:
timesns.append(round((f*dcdfreq*2e-6)+0.005,3))
print(timesns)
return(timesns)
loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/DMlipid/'
names=['05lp0','05lp40','05lp253']
#psimage.save('myImage.png')
# fig, axs =pl
fig, axes = plt.subplots(1, 4, figsize=(10, 10))
fig.set_figheight(7)
fig.set_figwidth(18)
axes.flatten()
row=0
tp=[0,40, 253]
sim='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/lipid24_05.dat'
data=get_time_dipole(sim,2500,1)
t=data[0]
theta=data[4]
axes[0].plot(t,theta)
axes[0].set_xlabel('Time(ns)',fontsize=Fntsize+4)
axes[0].set_ylabel('$\Theta_{z}$ (Degree)',fontsize=Fntsize+4)
timetpns=conv_frame_ns(tp,2500)
print(timetpns)
alltimens=timetpns
fieldname=['0.2V/nm','0.5V/nm']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)
for row, name in enumerate(names):
row +=1
psimage=Image.open(f'{loc}{name}')
axes[row].imshow(psimage)
axes[row].set_title(f'{alltimens[row-1]}ns',fontsize=Fntsize+5)
axes[row].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[row].xaxis.set_visible(False)
axes[row].tick_params(left=False, labelleft=False)
plt.setp(axes[row].spines.values(), visible=False)
if row==1:
axes[row].arrow(start_point[0], start_point[1],
end_point[0] - start_point[0], end_point[1] - start_point[1],
color='black', width=6, head_width=20, head_length=10)
axes[row].set_ylabel(f'$E_{z}$={fieldname[row]}',fontsize=20)
fig.suptitle('DW: 0.5V/nm',fontsize=30)
plt.tight_layout(pad=0.5)
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/lipidDM.png',dpi=300)
[0.005, 0.205, 1.27] [0.005, 0.205, 1.27]
from PIL import Image
def conv_frame_ns(frames,dcdfreq):
timesns=[]
for f in frames:
timesns.append(round((f*dcdfreq*2e-6)+0.005,3))
print(timesns)
return(timesns)
loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/Effectpepion/'
names=['DW05','D8CV02E','1DI02E']
#psimage.save('myImage.png')
# fig, axs =pl
fig, axes = plt.subplots(3, 4, figsize=(10, 10))
fig.set_figheight(18)
fig.set_figwidth(18)
#axes.flatten()
row=0
tp=[0,40, 253]
timetpns=conv_frame_ns(tp,2500)
print(timetpns)
alltimens=timetpns
fieldname=['0.5V/nm','0.2V/nm','0.2V/nm','0.2V/nm']
simname=['DW','D8CV', '1DI']
z='z'
start_point = (-1, 356)
end_point = (-1, 100)
times=[0,99,230,230 ]
alltimens=conv_frame_ns(times,2500)
for row, name in enumerate(names):
for col, frame in enumerate([0,99,230,230]):
if frame==230 and col==3:
psimage=Image.open(f'{loc}{name}_{frame}top')
print("hola")
else:
psimage=Image.open(f'{loc}{name}_{frame}')
axes[row,col].imshow(psimage)
if col==0 and row==0:
axes[row,col].set_title(f'{simname[row]}: {alltimens[col]}ns',fontsize=Fntsize+5)
elif col==0:
axes[row,col].set_title(f'{simname[row]}', fontsize=Fntsize+5)
elif row==0:
axes[row,col].set_title(f'{alltimens[col]}ns',fontsize=Fntsize+5)
axes[row,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[row,col].xaxis.set_visible(False)
axes[row,col].tick_params(left=False, labelleft=False)
plt.setp(axes[row,col].spines.values(), visible=False)
if col==0:
axes[row,col].arrow(start_point[0], start_point[1],
end_point[0] - start_point[0], end_point[1] - start_point[1],
color='black', width=6, head_width=20, head_length=10)
axes[row,col].set_ylabel(f'$E_{z}$={fieldname[row]}',fontsize=20)
plt.tight_layout(pad=0.5)
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/effecpepion.png',dpi=300)
[0.005, 0.205, 1.27] [0.005, 0.205, 1.27] [0.005, 0.5, 1.155, 1.155] hola hola hola
##########
## D8CV ##
##########
plt.rc('font', size=18)
fig, ax = plt.subplots(2, 2,figsize=(25, 15))
# fig = plt.figure()
# fig.set_figheight(14)
# fig.set_figwidth(20)
step=101
odr=9
ypt=0.9
Fntsize=30
lw=3
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
#simulation parameter
dif_con_sim=[4,2,5,3] #[4,1,6,3]
name_sim='0D8CV'
#Pmax=compute_max_normalzD8C(range(1,6),name_sim)
#Smoothing parameter
#plot parameters
transparency=0.3 #transparency of the data without simuation
#lw=3 #line with of the plot
labels=['0.0V/nm']
titles=["$E_{z}$", "-$E_{z}$"]
colors=['black']
lines=['solid']
#ax3 = ax.twinx()
#plot D8CV 0.0V/nm
for i,sim in enumerate([f'{location}{name_sim}_ZDMW.dat',f'{location}{name_sim}.dat']):
s=0
lw=2
x=np.array(get_time_dipole(sim,2500,1)[0])
y=np.array(get_time_dipole(sim,2500,1)[1])
theta=get_time_dipole(sim,2500,1)[4]
#suavizado
theta_smt=savgol_filter(theta,step,odr)
y_smt=savgol_filter(y,step,odr)
# ynorm_smt=y_smt/Pmax[i]
# ynorm=y/Pmax[i]
ynorm_smt=y_smt
ynorm=y
if i==0 :
ax[0,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[s], label=f'{labels[s]}',linewidth=lw)
ax[0,0].plot(x,ynorm,color=colors[s],linestyle=lines[s], alpha=0.3)
ax[0,1].plot(x,theta_smt,color=colors[s],label=f'{labels[s]}',linestyle=lines[s],linewidth=lw)
#axes[2,2].legend(bbox_to_anchor=(2.5,-0.2),fontsize=Fntsize+5,ncol=9)
#plot Pz and theta of DMPC
else:
ax[1,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[s],linewidth=lw,label=f'{labels[s]}')
ax[1,0].plot(x,ynorm,color=colors[s],linestyle=lines[s], alpha=0.3)
ax[1,1].plot(x,theta_smt,color=colors[s],linewidth=lw,linestyle=lines[s],label=f'{labels[s]}')
labels=[['0.1V/nm', '0.2V/nm'],['-0.1V/nm', '-0.2V/nm']]
colors=cm.rainbow(np.linspace(0, 1, 4))
lw=3
fields=[1,2]
tp=[99,68]
for row, namesimsufix in enumerate(['E','Er']):
for s , E in enumerate(fields):
#print("s:",s)
if E==1:
simulations=[f'{location}{name_sim}0{E}{namesimsufix}_ZDMW.dat',f'{location}{name_sim}0{E}{namesimsufix}_DMPC.dat']
else:
simulations=[f'{location}{name_sim}0{E}{namesimsufix}_ZDMW.dat',f'{location}{name_sim}0{E}{namesimsufix}.dat']
print(simulations)
lines=['solid','dotted']
a=0
for i in range(0,len(simulations)):
sim=simulations[i] #0.1V/nm
x=np.array(get_time_dipole(sim,2500,1)[0])
y=np.array(get_time_dipole(sim,2500,1)[1])
theta=get_time_dipole(sim,2500,1)[4]
#suavizado
theta_smt=savgol_filter(theta,step,odr)
y_smt=savgol_filter(y,step,odr)
# ynorm_smt=y_smt/Pmax[i]
# ynorm=y/Pmax[i]
ynorm_smt=y_smt
ynorm=y
#plot Pz and theta of water
if i==0 :
ax[0,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[row], label=f'{labels[row][s]}',linewidth=lw)
ax[0,0].plot(x,ynorm,color=colors[s],linestyle=lines[row], alpha=0.3)
ax[0,1].plot(x,theta_smt,color=colors[s],label=f'{labels[row][s]}',linestyle=lines[row],linewidth=lw)
#axes[2,2].legend(bbox_to_anchor=(2.5,-0.2),fontsize=Fntsize+5,ncol=9)
#plot Pz and theta of DMPC
else:
ax[1,0].plot(x,ynorm_smt,color=colors[s],linestyle=lines[row],linewidth=lw,label=f'{labels[row][s]}')
ax[1,0].plot(x,ynorm,color=colors[s],linestyle=lines[row], alpha=0.3)
ax[1,1].plot(x,theta_smt,color=colors[s],linewidth=lw,linestyle=lines[row],label=f'{labels[row][s]}')
if s==1 and i==1:
ax[0,0].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
ax[0,1].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
ax[1,0].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
ax[1,1].axvline(x=x[tp[row]], color=colors[s], linestyle=lines[row])
# ax[row,s].set_xlim(0,4)
ax[row,1].set_ylim(0,180)
#plot dipole:
from PIL import Image
def conv_frame_ns(frames,dcdfreq):
timesns=[]
for f in frames:
timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
return(timesns)
loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/dipoleMomentumDC/'
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
ax[0,0].set_ylabel('$P_{zw} (D)$',fontsize=Fntsize+5)
ax[1,0].set_ylabel('$P_{zDMPC} (D)$',fontsize=Fntsize+5)
#ax[2,0].set_ylabel('$P_{zall-pep} (D)$',fontsize=Fntsize+5)
ax[0,1].set_ylabel('$\Theta_{zw}$ (Degree)',fontsize=Fntsize+5)
ax[1,1].set_ylabel('$\Theta_{zDMPC}$ (Degree)',fontsize=Fntsize+5)
#ax[2,1].set_ylabel('$\Theta_{zall-pep}$ (Degree)',fontsize=Fntsize+5)
ax[1,0].set_xlim(0,6)
ax[1,1].set_xlim(0,6)
ax[1,0].set_xlabel('Time(ns)',fontsize=Fntsize+5)
ax[1,1].set_xlabel('Time(ns)',fontsize=Fntsize+5)
ax[1,0].legend(bbox_to_anchor=(2.1, -0.19),fontsize=Fntsize,ncol=6)
#plt.tight_layout(pad=3)
plt.show()
plt.show()
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/P_W_dmcp_p_D8CV.png',dpi=300)
chainpep={"DCV": ['A'], "D3CV": ['A','B','C'], "D8CV": ['A','B','C','D','E','F','G','H'],"DCH": ['A'], "D3CH": ['A','B','C']}
fig, axes = plt.subplots(4, 4, figsize=(10, 10))
fig.set_figheight(35)
fig.set_figwidth(40)
names_snap=['0D8CV01E','0D8CV01Er','0D8CV02E','0D8CV02Er']
titles=['D8CV: 0.1V/nm','D8CV: -0.1V/nm','D8CV: 0.2V/nm','D8CV: -0.2V/nm' ]
tppnD8CV=[99,68]
colors=['yellow','gray', 'orange','red', 'magenta', 'pink', 'green','indigo', 'black']
for col,name in enumerate(names_snap):
psimage=Image.open(f'{loc}{name}')
axes[0,col].imshow(psimage)
axes[0,col].set_title(f'{titles[col]}',fontsize=Fntsize+5)
axes[0,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[0,col].xaxis.set_visible(False)
axes[0,col].tick_params(left=False, labelleft=False)
plt.setp(axes[0,col].spines.values(), visible=False)
#axes[row,col].yaxis.set_visible(False)
sufixs=chainpep['D8CV']
for c,sufix in enumerate(sufixs):
sim=f'{location}{name}_p{sufix}.dat'
data=get_time_dipole(sim,2500,1)
t=data[0]
P=data[1]
theta=data[4]
axes[2,col].plot(t,P,label=f'Pep_{sufix}',color=colors[c])
axes[3,col].plot(t,theta,label=f'Pep{sufix}',color=colors[c])
names_snap=['0D8CV01Etop','0D8CV01Ertop','0D8CV02Etop','0D8CV02Ertop']
for col,name in enumerate(names_snap):
psimage=Image.open(f'{loc}{name}')
axes[1,col].imshow(psimage)
#axes[1,col].set_title(f'{titles[col]}',fontsize=Fntsize+5)
axes[1,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[1,col].xaxis.set_visible(False)
axes[1,col].tick_params(left=False, labelleft=False)
plt.setp(axes[1,col].spines.values(), visible=False)
axes[2,0].set_ylabel('$P_{zp} (D)$',fontsize=Fntsize+5)
axes[3,0].set_ylabel('$\Theta_{zp}$ (Degree)',fontsize=Fntsize+5)
for i in range(0,4):
axes[2,i].set_xlabel('Time(ns)',fontsize=Fntsize+5)
axes[3,2].legend(bbox_to_anchor=(2.4,-0.2),fontsize=Fntsize+5,ncol=9)
plt.show()
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/P_pnD8CV.png',dpi=300)
['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01E_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01E_DMPC.dat'] ['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02E_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02E.dat'] ['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01Er_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV01Er_DMPC.dat'] ['/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02Er_ZDMW.dat', '/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/0D8CV02Er.dat']
from PIL import Image
def conv_frame_ns(frames,dcdfreq):
timesns=[]
for f in frames:
timesns.append(round((f*dcdfreq*2e-6)+0.005,2))
return(timesns)
plt.rc('font', size=15)
loc='/media/thalia/New Volume/SIMULACIONES TDG/images TDG/dipoleMomentumDC/'
location='/media/thalia/New Volume/SIMULACIONES TDG/Electroporation_TDG/dipole momento/'
names_snap=['DCV','D3CV','D8CV','DCH','D3CH']
chainpep={"DCV": ['A'], "D3CV": ['A','B','C'], "D8CV": ['A','B','C','D','E','F','G','H'],"DCH": ['A'], "D3CH": ['A','B','C']}
colors=['yellow','gray', 'orange','red', 'magenta', 'pink', 'green','indigo' ]
fig, axes = plt.subplots(3, 5, figsize=(10, 10))
fig.set_figheight(20)
fig.set_figwidth(35)
for col,name in enumerate(names_snap):
psimage=Image.open(f'{loc}{name}')
axes[0,col].imshow(psimage)
axes[0,col].set_title(f'{name}',fontsize=Fntsize+5)
axes[0,col].tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False)
axes[0,col].xaxis.set_visible(False)
axes[0,col].tick_params(left=False, labelleft=False)
plt.setp(axes[0,col].spines.values(), visible=False)
#axes[row,col].yaxis.set_visible(False)
sufixs=chainpep[name]
for c,sufix in enumerate(sufixs):
sim=f'{location}{name}_p{sufix}.dat'
data=get_time_dipole(sim,2500,1)
t=data[0]
P=data[1]
theta=data[4]
axes[1,col].plot(t,P,label=f'Pep_{sufix}',color=colors[c])
axes[2,col].plot(t,theta,label=f'Pep_{sufix}',color=colors[c])
axes[1,0].set_ylabel('$P_{zp} (D)$',fontsize=Fntsize+5)
axes[2,0].set_ylabel('$\Theta_{zp}$ (Degree)',fontsize=Fntsize+5)
axes[2,2].set_xlabel('Time(ns)',fontsize=Fntsize+5)
axes[2,2].legend(bbox_to_anchor=(4,-0.2),fontsize=Fntsize+5,ncol=8)
plt.show()
# plt.tight_layout(pad=0.5)
fig.savefig('/media/thalia/New Volume/SIMULACIONES TDG/images TDG/pepDipoleDC.png',dpi=300)